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Colored noise and memory effects on formal spiking neuron models 
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Simplified neuronal models capture the essence of the electrical activity of a generic neuron, besides 
being more interesting from the computational point of view when compared to higher dimensional 
models such as the Hodgkin-Huxley one. In this work, we propose a generalized resonate-and-fire 
model described by a generalized Langevin equation that takes into account memory effects and 
colored noise. We perform a comprehensive numerical analysis to study the dynamics and the point 
process statistics of the proposed model, highlighting interesting new features like: i) non-monotonic 
behavior (emergence of peak structures, enhanced by the choice of colored noise characteristic time- 
scale) of the coefficient of variation (CV) as a function of memory characteristic time-scale, ii) colored 
noise-induced shift in the CV, and iii) emergence and suppression of multimodality in the interspike 
interval (ISI) distribution due to memory-induced subthreshold oscillations. Moreover, in the noise- 
induced spike regime, we study how memory and colored noise affects the coherence resonance (CR) 
phenomenon. We found that for sufficiently long memory, CR is not only suppressed, but also the 
minimum of the CV x noise intensity curve that characterizes the presence of CR may be replaced 
by a maximum. The aforementioned features allow to interpret the interplay between memory and 
colored noise as an effective control mechanism to neuronal variability. Since both variability and 
non-trivial temporal patterns in the ISI distribution are ubiquitous in biological cells, we hope the 
present model can be useful in modeling real aspects of neurons. 

PACS numbers: 87.19.11, 05.40.-a, 87.10.Mn, 87.19.1c 


I. INTRODUCTION 

The current consensus is that neurons are the fundamental units where information is processed in the nervous 
system. The celebrated Hodgkin-Huxley (HH) model mathematically translates into a four dimensional ordinary 
differential equation the detailed biophysical features of a neuron and nicely explains the generation of an action 
potential. Besides its relatively high dimensionality, another downside of the HH model rests on the fact that there 
are approximately twenty parameters to be determined. Consequently, a clear and intuitive understanding of the 
fundamentals of neuronal dynamics becomes difficult and one is often restricted to computer simulations. Also, this 
high dimensionality implies high computational costs when one aims at simulating a large population of neurons. 

If an accurate low-dimensional description of neuron activity could be found, it would certainly have advantages 
over high dimensional approaches, both in terms of simulation and understanding of the fundamental principles. In 
addition, the understanding of neuronal mechanisms resorting to simple models could be useful in the engineering 
of artificial neural devices designed to reproduce a given real biological feature. In that sense, several alternative 
approaches have been created along the years, like the spike-response class of models. Two notable examples in 
this class are the integrate-and-fire [I| and the resonate-and-fire models 0. These kinds of models are typically 
characterized by a differential equation whose solution describes the time evolution of the subthreshold membrane 
potential and an ad hoc rule according to which when the potential reaches the threshold (i) a spike is produced and 
(ii) the potential is reset to a certain value. Additionally, a refractory period, that is natural in the HH model, can 
be postulated in the such models. 

The precursor of integrate-and-fire model was introduced by the French physiologist Louis Lapicque in 1907, long 
before the introduction of the HH model Q. In 1965, Richard Stein introduced the so-called leaky integrate-and- 
fire (LIF) model which became very ptmular from the 1990s, when it started to be used (along side with its 
generalizations) in neural networks studies d, @ . Essentially, the integrate-and-fire model is constructed based on 
analogies with electrical circuits. We can picture the standard LIF model as a parallel resistive-capacitive (RC) 
circuit subjected to a deterministic external input. In addition to this deterministic input, there are several sources 
of fluctuation that influence the dynamics of the neuron potential. We can name at least three of these sources that 
stand out. The opening and closing of the ion channels behave stochastically Q, the release of neurotransmitters by 
chemical synapses is a Poissonian process, and the synaptic input from other neurons (of the order of 10^ synaptic 
junctions per neuron Q) can also be well modeled as a stochastic process Q. Stochastic versions of neuron models 


* leandro.silva@ufabc.edu.br 
^ rafael.vilela@ufabc.edu.br 



2 


have been strongly motivated by in vivo studies, since in this case the random time intervals of input synaptic impulses 
have to be taken into account. A prominent example in this class is the diffusive (or stochastic) leaky integrate-and-fire 
(SLIF). The main characteristic of the SLIF model is the explicit description of the randomness of synaptic inputs 
through the presence of a noise term in the equation of motion of the membrane potential v. 

Far from being just a mathematical device to dictate how the environment manifests itself in the effective dynamics 
of the system, the presence of noise in neural dynamics has proven fundamental in neurophysiological processes Q. 
Experimental evidence supports the idea that noise is intrinsically necessary to mechanisms responsible for signal 
detection, decision making, perception and short term memory [I EMI- Therefore, detailed studies connecting 
stochastic processes and neuronal dynamics at all levels (i.e. single neurons, networks and continuum field approaches) 
are necessary for a better understanding of the complexity of the nervous system. 

The subthreshold dynamics of the SLIF model is phenomenologically described by the following stochastic differ¬ 
ential equation, already written in the dimensionless form [l| : 

+v(t) = n + V^r]{t) , (1.1) 

where v{t) is the membrane potential at time t, r corresponds to the membrane time scale and /x is a constant average 
input if synaptic inputs are considered homogeneous. The resistor plays the role of the Ohmic leakage impedance of 
the neural membrane and the capacitor models the lipid bilayer that forms the neuronal membrane. In Eq. dm, 
r]{t) is to be interpreted as a Gaussian white noise term (i.e. (??(t)) = 0 and {'il{t)r]{t')) = S{t — t')) and a is the noise 
intensity. 

Eq. (11.11) corresponds to the simplest possible formulation of the SLIF model and may lose sight of some relevant 
biological features. In view of this, several generalizations have been proposed in the literature, many of them 
supported by experimental evidence. For example, Stevens and Zador suggest that the variability of output spike 
times of cortical neurons in the awake brain could be accounted for the existence of correlations in the arrival times of 
input at different synapses [l^. In Ref. (H), Salinas and Sejnowski highlight the necessity to take into account the 
impact of temporal correlations. Their motivation comes from experimental studies showing that cortical neurons can 
exhibit an input correlation time of the same order of magnitude as the mean interspike interval of the correspondent 
response. A practical way to implement a finite correlation time is to promote the term r]{t) in Eq. (|1.I[) to a colored 
noise term. 

Many other studies have been developed in order to achieve a more realistic (yet effective and simplified) approach 
to neuronal dynamics (T5l-[l8l|. and we would like to go a step forward. In a broader perspective, we can draw an 
analogy between neuronal dynamics and the nonequilibrium statistical theory of open systems. We can interpret 
a given neuron (or even a given neuronal tissue) in whose dynamics we are interested as the “system”. Since this 
system interacts with its surrounding medium, we have to account for the environment feedback into the system 
in order to obtain an accurate description of the system dynamics. In the early times of the development of the 
theory of open systems, the usual approach used to describe the effective dynamics of these systems was to model 
the problem resorting to a phenomenological Langevin equation. The original Langevin equation is characterized by 
Gaussian white noise and time-local dissipative terms. The main limitation of this approach is that the usual Brownian 
motion implicitly assumes that the environment interacts instantaneously with the system. This assumption is valid 
only in very restricted and particular limit cases. If a more realistic description of a given open system is required, 
it is necessary to take into account that information exchange (e.g. linear momentum exchange between massive 
particles, ionic current exchange between neurons and so on) takes place during finite time intervals, which must lead 
to finite memory effects and colored noise. In view of this, some phenomenological [l9j| and microphysical [20l - [^ 
approaches have been developed in order to elucidate how a more realistic effective equation of motion emerges when 
the environment degrees of freedom are coarse-grained. From this more fundamental perspective, a generalized theory 
of Brownian motion arose, which can be described by a generalized Langevin equation (GLE). In such an equation, 
nonlocal dissipative terms and correlated noise reflect finite-time effects. 

In this paper, we argue that it would be a step toward a more realistic modeling of effective neuronal dynamics 
to rewrite the equation that describes the membrane potential using the generalized Brownian motion framework. 
We thus propose to enrich intrinsically phenomenological neuron models to reflect a crucial aspect that characterizes 
first-principles models: the concomitant presence of memory effects and colored noise. This mathematical framework 
ensures causality and is particularly a requisite when the processes that take place in the system and environment 
have comparable characteristic time scales. We thus introduce and characterize a first-passage time model whose 
subthrehold dynamics obeys the Langevin equation of generalized Brownian motion. 

From the biophysical point of view, it is empirically known that the distribution of the interspike intervals (ISI) 
of real neurons (e.g. pallidal an d gang lion cells) can exhibit non-trivial patterns, like a bimodal or multimodal 
structure depending on the input (241 - 1^ . Neither the standard nor the stochastic LIF model is able to produce such 
a multimodal distribution. Even though the original resonate-and-fire model is capable of providing multimodality 
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under some conditions (28l - [30l| . it does not have an intrinsic dynamical mechanism to control the emergence and the 
shape of multimodality. We show that the memory term introduced by our model plays the role of a variability control 
mechanism. Adjusting the memory parameter, which reflects a property of the input, one enables smooth transitions 
between unimodal and multimodal regimes. 

The paper is organized as follows: in Sec.|TT]we introduce our model based on a first passage time (FPT) problem 
for the generalized Brownian motion, and review some analytical results on the subthreshold dynamics to guide the 
understanding of the numerical analysis, presented in Sec IIIII The results and discussions are the subject of Sec IIVI 
In Sec. E the outlook of this paper is given. 


II. A GENERALIZED RESONATE-AND-FIRE MODEL 

The SLIF model is an useful tool to mimic the potential dynamics of the so-called nonresonant cells. The response 
of this kind of cell to an external excitatory stimulus is characterized by a monotonic relaxation toward its resting 
state. This behavior becomes even clearer if we draw an analogy between Eq. dm and the equation of motion of 
an overdamped harmonic oscillator. In recent decades, a substantial amount of experimental data has shown that 
neurons can also exhibit damped subthreshold oscillations at a given frequency |3ll - l35l |. From a dynamical viewpoint, 
this occurs when the neurons operate close to an Andronov-Hopf bifurcation (36l | . The specific case in which there is no 
bistability and the rest state is the global attractor corresponds to a supercritical Andronov-Hopf bifurcation. Several 
models that aim at accounting for the aforementioned subthreshold oscillations have been constructed [H, [s!], [13,113 ■ 
Our starting point is an archetypal second-order model for the subthreshold dynamics, written in dimensionless 
form as: 


^ ^ , ( 2 . 1 ) 

where 77 (f) is again taken as white Gaussian noise. In terms of the analogy between neuronal dynamics and electrical 
circuits, the second-order derivative in the left hand side of Eq. (HID emerges when an inductor is plugged in parallel 
with both R and C in the LIF circuit. From the physiological point of view, the inductor emulates the presence 
of long-term processes, like calcium currents [3^ |4^. An important feature of the model given by Eq. (j2.1l) is that 
for 7 < 2u!^/X (underdamped case) the typical behavior of a resonant neuron is recovered. Analogously, 7 > 2a;'\/A 
corresponds to the overdamped case. In particular, in the limit case when A = 0 (usually referred to as the overdamped 
limit) the SLIF neuron that mimics the behavior of a nonresonant cell is recovered. Therefore, Eq. (12.11) can be thought 
of as a generalization of the SLIE model. Here we are interested in the case A 0. 

A more general approach corresponds to rewriting Eq. (EB using the framework of the generalized Brownian motion 
which, after rescaling the parameters in terms of A, reads: 


d?v{t) 

df^ 



dt'K{t - f) 


dv{t') 

dt' 


+ g'ivit))=m, 


( 2 . 2 ) 


where to is the time instant when the reset rule was last performed, K{t — t') is the memory kernel whose integral 
is normalized to the unit and whose characteristic time scale is denoted by 1/F, and ^ is a Gaussian colored noise 
with characteristic time scale 1/F^. The presence of colored noise is particularly relevant in situations where the 
characteristic time scales of fluctuations are comparable to the system time scale. The deterministic parameter g in 
Eq. (12.11) was absorbed in g'(v), which is a linear or nonlinear function of the membrane potential. The prime in g'{v) 
means derivative w.r.t v. If the neuron worked strictly as a RLG circuit and if the only source of fluctuation were 
of thermal origin, then a fluctuation-dissipation relation (EDR) would be expected to hold exactly (Nyquist noise): 
(^(t)C(t')) = ~ t')- Here, as mentioned above, we are motivated by other sources of noise as well, and do not 

impose beforehand such a perfect matching between noise and memory time scales. 

Equation (I2.2|l is the central equation of this paper. Once the memory kernel is defined (see Sec IIIII for a convenient 
choice), it describes the subthreshold evolution of the membrane potential. The establishment of a fire-and-reset rule 
completes the model. Here we adopt the rule: if v{t) = Wth then (i) a spike is considered to have occurred at time t and 
(ii) V ^ Vr and v' 0. Generally speaking, the reset rule reflects the biophysical fact that subsequent to a spike is 
the return of all the neuron state variables (e.g. in the Hodgkin-Huxley model) to a small region in phase space. The 
reset rule corresponds to the simplification that replaces that small region by a single point. Since the equation for 
the subthreshold dynamics, Eq. (12.21) . is of second-order for the variable v, it can be rewritten as a (two-dimensional) 
first-order equation for the variables v and dv/dt. These are the state variables of our model and, accordingly, must 
be both reset. The specific choice of the value zero for the reset of dv/dt does not have any special meaning and 
implies no loss of generality. 
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We can recover the standard first order SLIF model from Eq. (|2.2I1 by applying two approximations. The first one 
is the Markovian approximation, whose validity tends to be assured when the time scale of the system T = v/v is such 
that T ^ 1/r (additionally, T ^ l/hj if the white noise limit is also required), and for large times, when t/T ^ 1. 
The nonlocal term can be then rewritten as [l^, |4l| 


7 f dt'K{t — t')v{t') 

J Iq 


7 / dt' 26{t — t')v{t') = ^v{t) 


(2.3) 


At this point, we recover the second-order model Eq. (EtD with A = 1, which will be referred to as the Markovian 
approximation and compared to our proposed model Eq. (lO) in Sec. IIVI Going a step forward and performing a 
second approximation by imposing A = 0, we recover the hrst order SLIF model. The term g'(v) can be interpreted 
in mechanical terms as the derivative of an effective potential. For a linear choice of g'{v), we can define g{v) = 
12 — /iu, with minimum given by Vc = ■ For any v(t = 0) < Uth < Vc, given a sufficiently long time interval, 

the probability of spike occurrence is 1, both in the deterministic and in the stochastic cases. This regime is usually 
known as mean-driven or tonic firing regime. On the other hand, if v(t = 0) < i;c < the probability of occurrence 
of spikes in the overdamped regime is non-null only in the stochastic case. The same occurs in the underdamped 
regime if the maximum value of the time-dependent amplitude of oscillation around Vc is smaller than vtt when the 
noise is switched off. One refers to these cases as the fluctuation- or noise-induced firing regime. 

If we consider a linear version of (12.21) . standard Laplace transform techniques can be used to assess the 
subthreshold dynamics of the model |42l |. This procedure will be useful as a check for the numerical code and also 
to assist the FPT analysis in the next sections, when the fire-and-reset rule is turned on. Imposing a linear g'{v) 
function defined as g'{v) = uflv — /r in Eq. (USD, we obtain: 


v{s) 


where it was assumed to 


. h(0)+ [s-f 7 A:(s)] v{0) g{s) ^{s) 

-\-s'yK{s)-\-s'^-\-sjK{s)-\-s'^-\-s'yK(s)-\-uj'^ ' 

0 and the definition of the Laplace transform operator 


(2.4) 


^ {/}(«) = /(s) 



dt exp (—st)/(t) 


was used. The formal solution for the time evolution of the membrane potential is then given by: 

u(t) = {h} (t) = (/5(t)-I- f dt'h{t — t')g{t')-\- 

Jo 

where we have again resorted to the convolution theorem and dehned 


dt'h{t — t')^{t') 


(2.5) 


( 2 . 6 ) 


I -\- s^K{s) -I- I 


as well as 


h{t) = L 


-1 


W;- 


-I- s^K{s) Lvfl 

Noting that the noise term ^ is of zero mean, we can obtain the average dynamics of the potential v(t)\ 

poo 

{v{t)) = (p{t) -\- / dt'hit — t')/r(t') . 
do 

For the case of a constant g, Eq. (12.91) reduces to 

1/s 


(vit)) = ifit) -j- fiL 


-1 


-\- sjK{s) -\- up- 


it). 


(2.7) 


( 2 . 8 ) 


(2.9) 


( 2 . 10 ) 


Once the form of the memory kernel K{t — t') is specihed, we can easily evaluate Eq. (I2.10p numerically or alge¬ 
braically. For example, if we assume a Gamma kernel (which has as a particular case the exponentially decaying kernel) 
[ 4 ^ |43| or an exponentially damped harmonic kernel , analytical solutions are possible [12, ll^l • Since they 
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involve very long expressions, we refrain from writing them explicitly here. The standard deviation (T„ = sj (n^) — (n)^ 
of the membrane potential can be straightforward obtained considering Eq. (12.1011 and by noting that 




dt'h{t — 


dt"h{t-t") / dt'h{t-t'){^{t')i{t")) . (2.11) 

Jo 


Like in the first moment case, the second moment given by Eq. (12.111) can also be analytically obtained if, additionally 
to the memory kernel, a suitable correlation function for the noise term is provided. 

We point out that Karmeshu and Kadambari (Ref. [3) proposed a different approach to introduce memory ef¬ 
fects and also obtained interesting results. They generalized the standard first-order SLIF model by introducing a 
distributed delay in the following way: 


T pt 

= —/3 / dt'K{t — t')v{t') + ^1 + V2a"^r]{t) , 
dt Jq 

where r]{t) has the properties of a Gaussian white noise. 


( 2 . 12 ) 


III. NUMERICAL APPROACH AND SPIKE TRAIN STATISTICS 

A. Numerical scheme 

One of the main sources of information about the dynamics of a single neuron or an ensemble of neurons is the 
ISI distribution. The problem of obtaining such a distribution is formally known as a first passage time (FPT) 
problem [H, |4^ . Despite the indisputable success of the standard formulations on obtaining closed or approximated 
expressions for the FPT in the most diverse systems, the corresponding tools are not suitable to describe systems 
whose relaxation time is larger than the typical FPT [^, . This is precisely the case of the present model and also 

other models that are able to produce subthreshold oscillations. 

In Refs. Verechtchaguina et al. considered a resonant-and-fire model subjected to colored noise and found 

a closed form solution for the first term of the mean FPT equation, which in turn is formulated in terms of a sum 
of integrals. Despite the breakthrough developments of Verechtchaguina et ah, to improve the reliability of the 
approximation more terms need to be calculated semi-analytically, and then the computational cost for obtaining the 
approximated result rapidly approaches the one for simulating the dynamical equation. This is no less the case of the 
model introduced in this paper which, besides colored noise, also considers the presence of memory. For this reason, 
here we choose to perform an extensive numerical analysis of the problem. 

There are several methods in the literature to address integro-differential equations like Eq. (1^ For 

the memory kernel and noise correlation we use here, it is enough to adopt the method described in Ref. |52| . The 
objective is to map the nonlocal Langevin equation into a system of local stochastic differential equations. For this 
purpose, we introduce the auxiliary variable W, defined by 

W{t)=- f dt'K{t-t')v'(t') . (3.1) 

Jto 

We consider the simplest case for the memory kernel, a decaying exponential function also known as the Ornstein- 
Uhlenbeck kernel: 


=Texp[-r\t-t'\] . (3.2) 

For the noise term, we impose that its time evolution obey the Ornstein-Uhlenbeck stochastic differential equation: 

= ^ (3-3) 

where (t) is a white Gaussian noise satisfying 

iVi (t)) = 0 , 

(3.4) 

We use the stationary solution of Eq. (13.3p which, together with Eq. (13.4p , yields the following correlation function 
for the noise: 


(C (0 C it')) = exp [-Pf \t - t'\] . 


(3.5) 
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Now we can use the definition Eq. m and the differential equation Eq. (13.31) in order to rewrite the integro-differential 
Eq. (12.21) in terms of the following system of first order local differential equations: 

V {t) =y (t) 

y (t) =M - + iW (t) + ^{t) 

W{t) =-rW{t)- K{0)y(t) 

(t) - Y^2r|cr| 775 (0 ■ (3-6) 

The numerical model is completed upon a specification of the fire-and-reset rule. Here it reads: if v(t) = uth then 
(i) a spike is considered to have occured at time t and (ii) u —>■ Ur, ?/ —>■ 0, and VE —>■ 0. 

We have used a standard fourth-order Runge-Kutta method with a time step size varying between At = 10“^ and 
At = 10“^ to obtain the numerical results shown in the following. This numerical scheme was found to be enough 
both for numerical stability and numerical precision. 


B. Standard measuring tools 

Neuronal dynamics can be characterized by several quantities [l|, Is^ - lssjl and here we define some of them that will 
be useful in the analysis we perform in the next sections. The spike time is defined as the time instant when the 
membrane potential reaches the threshold value. Given a time window = tf — tg and a trial (say, the Hh-trial), 
a neuron produces a finite set of spikes, whose temporal distribution can be accounted for by 

Ni 

y^it) = ^ (3-7) 

with t\ > to and The distribution yi{t) is called spike train, and the sum runs discretely over each spike 

time t*, from j = 1 to j = iV^. A consistency relation between yiit) and Ni can be written. The number of spikes in 
a given trial can be evaluated as 


ftf 

Ni= dtyi{t) . 
Jto 


(3.8) 


In n trials, the total number of spikes is defined as N = X)r=i 


We define the firing rate as an average over time: 


N 

nAt\Y 


(3.9) 


The time interval between two successive spikes ( Tj = — tj) is called the interspike interval (ISI). This 

definition allows us to concatenate the ISPs of all n trials and then define the following quantities: the average ISI 
(denoted by (T)), 


(T) 


1 

N-n 


n Ni — 1 

EEr‘. 


i=l j=l 


and the coefficient of variation (CV), denoted by c„ in the equations: 


Cy 


(T) 


(3.10) 


(3.11) 


where (T^) is defined analogously to Ea. (l3.10[) . The CV is of course a measure of the regularity of the spike train. 
From Eq. (I3.11|) . one can infer that a time uniform spike train leads to = 0, while a Poissonian spike train is 
characterized by Ct, = 1. 
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Markovian case 



c =0.43 

V 
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N = 9 
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V 
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Non-Markovian, F = = 0.3 



t 


Figure 1. Sample trials of the potential v{t) obtained from Eg. 113.611 in the tonic spiking regime using = 0.1. 


IV. RESULTS AND DISCUSSIONS 

In the following numerical simulations, we have considered statistics based on at least n = 10^ trials where each 
trial is characterized by an observational time window Atw = 10^. These choices were made in order to obtain a 
sample space composed by at least 10® spikes. In all cases, we set a threshold value Uth = 0.1 and a reset value 
Vr = —0.05. We distinguish between two characteristic regimes: in subsection IIV Al we analyze the neuron dynamics 
when > Uth (tonic spiking regime, /i = 0.2), and in subsection IIV Bl we analyze the situation where /i/w^ < Uth, 
the so-called noise-induced spiking regime = 0.08). We considered in all cases 7 = 5 and w = 1, corresponding to 
the overdamped case. Therefore, when present, the subthreshold oscillations are due to the presence of memory. 


A. Tonic spiking regime 

We start by illustrating the dynamics provided by Eq. (13.61) showing in Fig. [T] some sample trials in a reduced time 
window for different memory parameter T. 

In Fig. [21 we test the limiting behavior of our numerical code. As expected, for large values of F and T^, we observe 
a good agreement between the non-Markovian dynamics and the Markovian approximation, given by Eq. with 
A = 1: on the left panel of Fig. [2] we show this agreement for the probability density function of the interspike intervals. 
On the right panel, we show the approach of the coefficient of variation to its counterpart in the Markovian limit for 
large values of F. 

In the inset plot of Fig. |21 we call attention to the non-monotonic behavior of the non-Markovian CV as a function 
of the memory parameter. In the long memory limit (F = F^ <C 1), the CV is near zero, what characterizes a very 
regular spike train. If we proceed increasing F, the value of CV also increases until a local maximum is reached 
(F = F^ Ri 0.5). After that, the value of CV starts to decrease until it reaches the Markovian limit (F = F^ ^ I). 
In Fig. [3] (left panel) we observe that increasing the noise intensity shifts the CV curves toward a higher level of 
variability, and also causes a softening effect on the shape of the CV peak. On the right panel of the same figure, we 
show the impact of noise on the mean firing rate r: increasing results in an appreciable decrease of the firing rate 
in the approximated range 0.I<F = F^<0.8, and in an increase of the firing rate for F = F^ > 0.8, approximately. 

The non-trivial behavior induced by memory described above, which is also present in the noise-induced spike 
regime (Sec.IV B), can be seen as an effective representation of some underlying mechanism that acts controlling the 
firing variability and rate of a given neuron. A number of such mechanisms have been reported in the literature. 
For instance, it has been suggested that GABAergic autaptic transmission (GAT) is responsible for the regulation of 

















































































Figure 2. (Color online) Convergence to the Markovian dynamics. Left: probability density function for the interspike interval 
obtained using Eg. (13.611 and its Markovian approximation. Right: the coefficient of variation (CV) as a function of F. The inset plot 
shows a detailed view of the CV for shorter values of F. Note the nonmonotonic behavior of CV as we increase F CMarkovian case: 




Figure 3. The behavior of the CV and the firing rate r as a function of F = F^ in the tonic spiking regime (p = 0.2). 


spike-timing precision in the so-called fast spiking interneurons, leading to the suppression of variability . In 

the case of our model, the memory is introduced to reflect a property of the input. This has a strong analogy with 
the experimental results reported in Refs, and [H^, which show that the overall activity level of the network that 
embeds a given neuron modifies the integrative properties and the performance of that single cell. 

It is important to note that it is the presence of memory rather than that of colored noise the determinant factor 
that induces the emergence of a peak in the c„x F curve. If we ease the restriction F = Fj (used only to reduce the 
number of free parameters) and consider Eq. (13.61) in the white noise limit (Fj —>■ oo), we obtain the result shown in 
the left panel of Fig. 2] the dotted line, that essentially coincides with the Fj = 100 curve (full line), still exhibits 
the local maximum of the CV as a function of F. Note that, in sharp contrast, no maximum occurs for the CV as 
a function of F^ when we consider the memoryless limit (F —^ oo) of the dynamical equation (approximated by the 
F = 100 curve on the right panel of Fig. n. Looking again at the left panel of Fig. 21 one notes that the value of 
Fj determines how smooth the peak structure is. The smaller the value of Fj is, the sharper is the peak, and even 
a double-peak structure can emerge in the Ci, x F curve (dashed line, Fj = 0.1). The bigger the value of Fj is (until 
it saturates in the white noise limit), the smoother is the shape of the peak. Also, increasing F^ implies a global 
upward shift in the neuron variability, situation that tends to be slightly violated only in the vicinity of the peak 
when Fj ^ 1, since the peak becomes increasingly smooth as Fj is increased. 

A joint analysis of the subthreshold evolution of Eq. (13.6p and the ISI histograms for a given time scale F“^ is 






































































Figure 4. Dissociated effects of memory and colored noise. Left: CV as a function of F for different values of F^. Right: 
Suppression of the peak in the memoryless limit, F 3> 1. In both panels, p = 0.2 and crj = 0.1. 


useful to shed light on the mechanism that leads to the emergence of the peak in the x F curve. In Fig. [5] we take 
a closer look at the case = 0.1 by showing a sequence of snapshots for fixed values of F = F^. In panel (a), we 
set F = F^ = 0.3, a value to the left of the maximum F = F^ Ri 0.5 (see Fig. El), and it can be noted that the ISI 
distribution exhibits a very sharp unimodal structure, which is characterized by a small value of CV, as shown in Fig.|3l 
These findings can be easily understood if one looks at the inset plot, which shows the mean membrane potential 
trajectory (Eq. 12.101) crossing the threshold value at a very early time instant due to memory-induced oscillations. 
The amplitude of these oscillations is so large that non-crossing trajectories become unlikely already after the first 
peak of the mean potential. When F = F^ = 0.5 (panel (b)), we can note a decrease in the height of the first ISI peak 
and the emergence of a second one. This configuration of the ISI distribution characterizes the maximum CV, since 
trajectories that are able to avoid the threshold crossing at the first oscillation do not possess a significantly dominant 
most probable interspike interval, since the emergent mode does not exhibit a sufficiently pronounced peak. Setting 
F = F^ = 0.8 (panel (c)), the second mode becomes well defined, since its height is now comparable to the one of 
the first peak. Thus, in this configuration we can infer a decrease in the CV value, what is corroborated by Fig. [31 
This decreasing behavior of the CV persists from this panel to the last one. Note that for larger values of F = F^, we 
have a continuous suppression of the bimodality, as can be seen in panels (e) and (f). However, even though both (a) 
and (f) ISI distributions exhibit a single-peaked structure, they yield different values of the CV: as the system moves 
toward the Markovian limit (F = F^ ^ 1) plateau exhibited in the right panel of Fig. [3J the CV decreases, but it 
saturates in a value larger than the one obtained when F = F^ <C 1. 

The same analysis can be performed to clarify the double-peak structure in the Ci, x F curve that is exhibited in the 
F^ = 0.1 case of Fig. [d] In Fig. [6] we show the ISI distribution for representatives values of F. From panel (a) to panel 
(b), we note the transition from a sharp, unimodal distribution to a distribution with two dominant modes. Panel 
(b) represents the maximum CV configuration for this particular case. In panel (c), we fix F = 0.7 (the minimum 
between the two peaks in Fig. HI) and a suppression of the first mode occurs, thus decreasing the CV value since only 
one dominant mode remains. The next panel corresponds to the second peak shown in Fig. 11, and from this panel 
to the last one, we note a decrease in the CV value until the memory less limit is reached and the CV assumes an 
approximately constant value. 

As a last remark in this subsection, we illustrate in Fig. |7| how the coefficient of variation behaves as a function 
of the noise intensity a^. The results are intuitive, since we can note a monotonic increase in the value of the CV 
as the noise intensity increases. The dotted red curves that are present in all panels correspond to the Markovian 
approximation (F = F^ —>• oo). As expected, for large values of F = F^ (panel (c), F^ = 100 case), the Markovian 
dynamics becomes a reliable approximation for the non-Markovian dynamics. 


B. Noise-induced spike regime and coherence resonance 

Now let us turn to the analysis of the noise-induced spike regime, where the external input /i and the parameter uj 
are defined such that pL/ta'^ < rith. 

In this regime, the qualitative behavior of both the CV and rate as functions of (F = F^) is similar to their 


































Figure 5. Emergence of bimodality in the ISI distribution. In all panels, /r = 0.2 and ctj = 0.1. 


counterparts in the tonic spiking regime (Fig. [TO] However, if we take a closer look at the ISI distribution in this case 
(Fig. ini where the = 0.1 case of Fig. [TUll is depicted), we can note the emergence of a multimodal structure that 
results from the memory induced subthreshold oscillations. Like in the previous subsection, for F = Fj <C 1 (long 
memory limit, panel (a)) and for F = Fj ^ 1 (Markovian limit, panel (f)), the ISI distribution is unimodal, and 
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Figure 6. Emergence of bimodality in the ISI distribution: analysis of the double-peak structure in the = 0.1 case of Fig. [4] 
In all panels: Fj = 0.1, p = 0.2 and = 0.1. 


for a critical value F = Fj Ri 0.1 (panel (b), log scale), we obtain a multimodal ISI distribution that maximizes the 
value of CV. It is important to note that we adopted the standard definition of the noise-induced spike regime, thus 
taken from the memoryless model perspective. Since memory is able to induce oscillatory behavior in the potential 


























































































































































































Figure 7. (Color online) The behavior of the CV as a function of the additive noise intensity in the tonic spiking regime. 


dynamics, it can eventually cause deterministic threshold crossing depending on the amplitude of the oscillations. For 
example, panels (a) and (b) show sufficiently small values of F = Fj such that memory provides a transition from 
noise-induced to tonic spike regime. From panel (c) onward, the spike regime becomes noise-induced in the standard 
sense. Also, from panel (b) to panel (e), we observe a decrease in the CV value until it reaches the Markovian limit 
value in panel (f). In this panel, we recover the typical unimodal ISI histogram of overdamped models. 

Now we discuss an interesting phenomenon known as coherence resonance (CR) [b^, that can be present not only 
in neuronal systems [63, but also in optical systems, plasma discharge dynamics and others [63 - 1^ . Like in 
stochastic resonance (SR) |66l. 1^. CR can be interpreted as a non-deleterious behavior induced by the presence of 
noise. But unlike SR, where setting an appropriate amount of noise can improve the ability of neurons detect external 
subthreshold signals (e.g. periodic inputs), the constructive role of noise in CR phenomenon is characterized by the 
emergence of order under the presence of a stochastic source only. Specifically in the neuronal dynamics case, a 
critical level of noise intensity can maximize the regularity (i.e. the coefficient of variation reaches a minimum) of 
a spike train. We now investigate how the introduced memory effects and colored noise affects the CR mechanism. 
For that purpose, in Fig. [10] we plot the CV as a function of the noise intensity. The dotted red line in plot (a) 
shows the result for the Markovian limit (F = Fj ^ 1) of Eq. (13.61) . An interesting result is that for sufficiently 
long memory (F = Fj = 0.1 case), CR is suppressed: the minimum in the Cy x curve is replaced by a maximum. 
Therefore, instead of the existence of a critical value for such that the ISI regularity is maximized, for long memory 
we have a critical value of the noise intensity such that the ISI regularity is minimized. Therefore, memory is able to 
not only suppress CR, but it can also induce an opposite situation, thus characterizing an “incoherence resonance” 
phenomena. This behavior is a manifestation of the deterministic threshold crossing induced by memory oscillations, 
as aforementioned. For the curves with F = Fj = 0.1 and 0.05, the oscillation amplitude is sufficiently large such 






























Figure 8. The behavior of the CV and the firing rate r as a function of the F = in the noise induced spike regime. 


that a transition to the tonic spike regime occurs, then completely changing the qualitative aspect of the c„ x 
curve. It is worth noting that the authors of Ref. [6lj| also describe an interesting coherence-incoherence maximization 
mechanism regulated by the noise level in the context of a SLIF that includes an absolute refractory period. 

In Fig. fTOT b'l. we call attention to the fact that, again, memory shows up as the fundamental ingredient that leads 
to the qualitative change in the c„ x cr^, since the characteristic time scale of the colored noise is not able to change 
the qualitative aspect of the curves, nor in the F^ <C 1 cases neither in F^ 1 cases. 


V. CONCLUSIONS AND FINAL REMARKS 


We have proposed a generalization of the stochastic resonate-and-fire model to the context of the generalized Brow¬ 
nian motion framework. We analyzed in detail the consequences of the concomitant presence of memory (distributed 
delay) and colored noise on the membrane potential dynamics. The model shows a remarkably rich dynamics, ex¬ 
hibiting a non-trivial behavior of the coefficient of variation as a function of the memory characteristic time-scale 1/F, 
thus suggesting that memory can be seen as an effective mechanism for generating and controlling neuron variability. 
Also, the present model is able to produce ISI distributions that exhibit multimodality, a feature whose temporal pat¬ 
tern can be controlled by adjusting the memory time-scale, then allowing smooth transitions between unimodal and 
multimodal distributions. This is a very desirable mechanism, since real neurons are able to exhibit such non-trivial 
patterns in the ISI distribution. Additionally, we have discussed the role played by colored noise in the membrane 
potential dynamics. We have shown that i) by adjusting the noise intensity it is possible to shift the x (F = F^) 
curve and to smooth its peak structure, and ii) by adjusting the noise time-scale I/Fj a double-peak structure can 
emerge in the Ct, x F curve. Finally, we have studied how memory and colored noise effects modifies the establishment 
of the coherence resonance phenomenon. We found that long memory is able to suppress the presence of CR, and 
for a given range of F and F^, a maximum in the Cy x is present instead of the usual minimization of CV that 
characterizes the coherence resonance. Given the variety of dynamical regimes provided by the proposed theoretical 
framework, we think that the present formulation could be useful to shed some new light on the modeling of dynamical 
aspects of real neurons. Also, the results presented here could be promptly extended to the context of other kind of 
noises, like the I//, for example. 
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Figure 9. Emergence of multimodality in the ISI distribution. In all panels: jj, = 0.08 and ctj = 0.1. 

















































































































Figure 10. (Color online) The behavior of the CV as a function of the noise intensity cr^ in the noise-induced spike regime: 
coherence resonance can be suppressed by memory effects. 
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